{
 "metadata": {
  "name": "",
  "signature": "sha256:31a58958351802834d98e0d3b9e17b7629fd18f7d1cc0d5af48ce95fe55a7eb4"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Spherical functions"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from __future__ import print_function, division, absolute_import\n",
      "\n",
      "import numpy as np\n",
      "import random\n",
      "import quaternion\n",
      "import numbapro as nb\n",
      "from numbapro import *\n",
      "import spherical_functions as sp\n",
      "\n",
      "ru = lambda : random.uniform(-1,1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Wigner 3-j symbols"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "j1,j2,j3,m1,m2,m3 = 1, 0, 1, 0, 0, (-0 - 0)\n",
      "sp.Wigner3j(j1,j2,j3,m1,m2,m3), sp.Wigner3j(j1,(j2+j3-m1)/2,(j2+j3+m1)/2,j3-j2,(j2-j3-m1)/2+m1+m2,(j2-j3+m1)/2-m1-m2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 2,
       "text": [
        "(-0.5773502691896257, 0.0)"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "j1,(j2+j3-m1)/2,(j2+j3+m1)/2,j3-j2,(j2-j3-m1)/2+m1+m2,(j2-j3+m1)/2-m1-m2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 3,
       "text": [
        "(1, 0.5, 0.5, 1, -0.5, -0.5)"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "j2,j3,m1,(j2-j3-m1)/2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "(0, 1, 0, -0.5)"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%timeit\n",
      "for j1 in range(8):\n",
      "    for j2 in range(8):\n",
      "        for j3 in range(8):\n",
      "            for m1 in range(-j1,j1+1):\n",
      "                for m2 in range(-j2,j2+1):\n",
      "                    sp.wigner3j(j1,j2,j3,m1,m2,-m1-m2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10 loops, best of 3: 22.5 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "i=0\n",
      "for j1 in range(8):\n",
      "    for j2 in range(8):\n",
      "        for j3 in range(8):\n",
      "            for m1 in range(-j1,j1+1):\n",
      "                for m2 in range(-j2,j2+1):\n",
      "                    i+=1\n",
      "i"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 3,
       "text": [
        "32768"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "22.5e-3/_"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "6.866455078125e-07"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Wigner $\\mathfrak{D}$ symbols"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "q = np.quaternion(ru(), ru(), ru(), ru())\n",
      "q = q/q.abs()\n",
      "q"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 2,
       "text": [
        "quaternion(0.542904926009884, -0.513641541910117, 0.265393949689424, -0.609091667326559)"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Using slow wrapper\n",
      "\n",
      "ells = range(16+1)\n",
      "evals = np.empty_like(ells, dtype=int)\n",
      "nanoseconds = np.empty_like(ells, dtype=float)\n",
      "for i,ell_max in enumerate(ells):\n",
      "    indices = np.array([[ell, mp, m] for ell in range(ell_max+1) for mp in range(-ell,ell+1) for m in range(-ell,ell+1)], dtype=int)\n",
      "    elements = np.empty((indices.shape[0],), dtype=complex)\n",
      "    result = %timeit -o sp.WignerD(q, indices)\n",
      "    evals[i] = len(indices)\n",
      "    nanoseconds[i] = 1e9*result.best/len(indices)\n",
      "    print(\"With ell_max={0}, and {1} evaluations, each D component averages {2:.0f} ns\".format(ell_max,evals[i],nanoseconds[i]))\n",
      "    sys.stdout.flush()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100000 loops, best of 3: 5.55 \u00b5s per loop\n",
        "With ell_max=0, and 1 evaluations, each D component averages 5555 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10000 loops, best of 3: 21.9 \u00b5s per loop\n",
        "With ell_max=1, and 10 evaluations, each D component averages 2194 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10000 loops, best of 3: 67.2 \u00b5s per loop\n",
        "With ell_max=2, and 35 evaluations, each D component averages 1920 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10000 loops, best of 3: 154 \u00b5s per loop\n",
        "With ell_max=3, and 84 evaluations, each D component averages 1833 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 295 \u00b5s per loop\n",
        "With ell_max=4, and 165 evaluations, each D component averages 1786 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 522 \u00b5s per loop\n",
        "With ell_max=5, and 286 evaluations, each D component averages 1825 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 825 \u00b5s per loop\n",
        "With ell_max=6, and 455 evaluations, each D component averages 1813 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 1.28 ms per loop\n",
        "With ell_max=7, and 680 evaluations, each D component averages 1877 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 1.77 ms per loop\n",
        "With ell_max=8, and 969 evaluations, each D component averages 1830 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 2.32 ms per loop\n",
        "With ell_max=9, and 1330 evaluations, each D component averages 1746 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 3.07 ms per loop\n",
        "With ell_max=10, and 1771 evaluations, each D component averages 1731 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 4.03 ms per loop\n",
        "With ell_max=11, and 2300 evaluations, each D component averages 1750 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 4.99 ms per loop\n",
        "With ell_max=12, and 2925 evaluations, each D component averages 1705 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 6.16 ms per loop\n",
        "With ell_max=13, and 3654 evaluations, each D component averages 1686 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 8.07 ms per loop\n",
        "With ell_max=14, and 4495 evaluations, each D component averages 1796 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 9.39 ms per loop\n",
        "With ell_max=15, and 5456 evaluations, each D component averages 1722 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 11.6 ms per loop\n",
        "With ell_max=16, and 6545 evaluations, each D component averages 1772 ns\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Directly using jitted inner function, also creating array\n",
      "\n",
      "ells = range(16+1)\n",
      "evals = np.empty_like(ells, dtype=int)\n",
      "nanoseconds = np.empty_like(ells, dtype=float)\n",
      "for i,ell_max in enumerate(ells):\n",
      "    indices = np.array([[ell, mp, m] for ell in range(ell_max+1) for mp in range(-ell,ell+1) for m in range(-ell,ell+1)], dtype=int)\n",
      "    result = %timeit -o sp._WignerD(q.a, q.b, indices, np.empty((indices.shape[0],), dtype=complex))\n",
      "    evals[i] = len(indices)\n",
      "    nanoseconds[i] = 1e9*result.best/len(indices)\n",
      "    print(\"With ell_max={0}, and {1} evaluations, each D component averages {2:.0f} ns\".format(ell_max,evals[i],nanoseconds[i]))\n",
      "    sys.stdout.flush()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000000 loops, best of 3: 1.28 \u00b5s per loop\n",
        "With ell_max=0, and 1 evaluations, each D component averages 1284 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100000 loops, best of 3: 3.06 \u00b5s per loop\n",
        "With ell_max=1, and 10 evaluations, each D component averages 306 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100000 loops, best of 3: 8.69 \u00b5s per loop\n",
        "With ell_max=2, and 35 evaluations, each D component averages 248 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10000 loops, best of 3: 20.3 \u00b5s per loop\n",
        "With ell_max=3, and 84 evaluations, each D component averages 242 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10000 loops, best of 3: 41.6 \u00b5s per loop\n",
        "With ell_max=4, and 165 evaluations, each D component averages 252 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10000 loops, best of 3: 74.3 \u00b5s per loop\n",
        "With ell_max=5, and 286 evaluations, each D component averages 260 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10000 loops, best of 3: 122 \u00b5s per loop\n",
        "With ell_max=6, and 455 evaluations, each D component averages 268 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10000 loops, best of 3: 186 \u00b5s per loop\n",
        "With ell_max=7, and 680 evaluations, each D component averages 274 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 271 \u00b5s per loop\n",
        "With ell_max=8, and 969 evaluations, each D component averages 280 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 381 \u00b5s per loop\n",
        "With ell_max=9, and 1330 evaluations, each D component averages 286 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 515 \u00b5s per loop\n",
        "With ell_max=10, and 1771 evaluations, each D component averages 291 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 680 \u00b5s per loop\n",
        "With ell_max=11, and 2300 evaluations, each D component averages 295 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 877 \u00b5s per loop\n",
        "With ell_max=12, and 2925 evaluations, each D component averages 300 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 1.12 ms per loop\n",
        "With ell_max=13, and 3654 evaluations, each D component averages 307 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 1.37 ms per loop\n",
        "With ell_max=14, and 4495 evaluations, each D component averages 305 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 1.65 ms per loop\n",
        "With ell_max=15, and 5456 evaluations, each D component averages 303 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 1.99 ms per loop\n",
        "With ell_max=16, and 6545 evaluations, each D component averages 304 ns\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# After re-configuring the calculation\n",
      "\n",
      "ells = range(16+1)\n",
      "evals = np.empty_like(ells, dtype=int)\n",
      "nanoseconds = np.empty_like(ells, dtype=float)\n",
      "for i,ell_max in enumerate(ells):\n",
      "    indices = np.array([[ell, mp, m] for ell in range(ell_max+1) for mp in range(-ell,ell+1) for m in range(-ell,ell+1)], dtype=int)\n",
      "    elements = np.empty((indices.shape[0],), dtype=complex)\n",
      "    result = %timeit -o sp._WignerD(q.a, q.b, indices, elements)\n",
      "    evals[i] = len(indices)\n",
      "    nanoseconds[i] = 1e9*result.best/len(indices)\n",
      "    print(\"With ell_max={0}, and {1} evaluations, each D component averages {2:.0f} ns\".format(ell_max,evals[i],nanoseconds[i]))\n",
      "    sys.stdout.flush()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000000 loops, best of 3: 558 ns per loop\n",
        "With ell_max=0, and 1 evaluations, each D component averages 558 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100000 loops, best of 3: 2.04 \u00b5s per loop\n",
        "With ell_max=1, and 10 evaluations, each D component averages 204 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100000 loops, best of 3: 7.18 \u00b5s per loop\n",
        "With ell_max=2, and 35 evaluations, each D component averages 205 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100000 loops, best of 3: 17.7 \u00b5s per loop\n",
        "With ell_max=3, and 84 evaluations, each D component averages 211 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10000 loops, best of 3: 35.2 \u00b5s per loop\n",
        "With ell_max=4, and 165 evaluations, each D component averages 213 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10000 loops, best of 3: 66.7 \u00b5s per loop\n",
        "With ell_max=5, and 286 evaluations, each D component averages 233 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10000 loops, best of 3: 109 \u00b5s per loop\n",
        "With ell_max=6, and 455 evaluations, each D component averages 239 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10000 loops, best of 3: 165 \u00b5s per loop\n",
        "With ell_max=7, and 680 evaluations, each D component averages 243 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 243 \u00b5s per loop\n",
        "With ell_max=8, and 969 evaluations, each D component averages 250 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 336 \u00b5s per loop\n",
        "With ell_max=9, and 1330 evaluations, each D component averages 252 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 456 \u00b5s per loop\n",
        "With ell_max=10, and 1771 evaluations, each D component averages 258 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 620 \u00b5s per loop\n",
        "With ell_max=11, and 2300 evaluations, each D component averages 270 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 792 \u00b5s per loop\n",
        "With ell_max=12, and 2925 evaluations, each D component averages 271 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 1.01 ms per loop\n",
        "With ell_max=13, and 3654 evaluations, each D component averages 277 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 1.25 ms per loop\n",
        "With ell_max=14, and 4495 evaluations, each D component averages 279 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 1.57 ms per loop\n",
        "With ell_max=15, and 5456 evaluations, each D component averages 287 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 1.87 ms per loop\n",
        "With ell_max=16, and 6545 evaluations, each D component averages 285 ns\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Directly using jitted inner function\n",
      "\n",
      "ells = range(16+1)\n",
      "evals = np.empty_like(ells, dtype=int)\n",
      "nanoseconds = np.empty_like(ells, dtype=float)\n",
      "for i,ell_max in enumerate(ells):\n",
      "    indices = np.array([[ell, mp, m] for ell in range(ell_max+1) for mp in range(-ell,ell+1) for m in range(-ell,ell+1)], dtype=int)\n",
      "    elements = np.empty((indices.shape[0],), dtype=complex)\n",
      "    result = %timeit -o sp._WignerD(q.a, q.b, indices, elements)\n",
      "    evals[i] = len(indices)\n",
      "    nanoseconds[i] = 1e9*result.best/len(indices)\n",
      "    print(\"With ell_max={0}, and {1} evaluations, each D component averages {2:.0f} ns\".format(ell_max,evals[i],nanoseconds[i]))\n",
      "    sys.stdout.flush()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000000 loops, best of 3: 615 ns per loop\n",
        "With ell_max=0, and 1 evaluations, each D component averages 615 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100000 loops, best of 3: 2.26 \u00b5s per loop\n",
        "With ell_max=1, and 10 evaluations, each D component averages 226 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100000 loops, best of 3: 7.87 \u00b5s per loop\n",
        "With ell_max=2, and 35 evaluations, each D component averages 225 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100000 loops, best of 3: 19.5 \u00b5s per loop\n",
        "With ell_max=3, and 84 evaluations, each D component averages 233 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10000 loops, best of 3: 41.5 \u00b5s per loop\n",
        "With ell_max=4, and 165 evaluations, each D component averages 251 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10000 loops, best of 3: 72.8 \u00b5s per loop\n",
        "With ell_max=5, and 286 evaluations, each D component averages 255 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10000 loops, best of 3: 120 \u00b5s per loop\n",
        "With ell_max=6, and 455 evaluations, each D component averages 264 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10000 loops, best of 3: 186 \u00b5s per loop\n",
        "With ell_max=7, and 680 evaluations, each D component averages 273 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 269 \u00b5s per loop\n",
        "With ell_max=8, and 969 evaluations, each D component averages 278 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 381 \u00b5s per loop\n",
        "With ell_max=9, and 1330 evaluations, each D component averages 287 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 520 \u00b5s per loop\n",
        "With ell_max=10, and 1771 evaluations, each D component averages 294 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 686 \u00b5s per loop\n",
        "With ell_max=11, and 2300 evaluations, each D component averages 298 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 876 \u00b5s per loop\n",
        "With ell_max=12, and 2925 evaluations, each D component averages 300 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 1.12 ms per loop\n",
        "With ell_max=13, and 3654 evaluations, each D component averages 307 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 1.37 ms per loop\n",
        "With ell_max=14, and 4495 evaluations, each D component averages 304 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 1.65 ms per loop\n",
        "With ell_max=15, and 5456 evaluations, each D component averages 302 ns\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 1.98 ms per loop\n",
        "With ell_max=16, and 6545 evaluations, each D component averages 303 ns\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Before wrapping key math parts in jit\n",
      "\n",
      "ells = range(1,16)\n",
      "evals = np.empty_like(ells, dtype=int)\n",
      "microseconds = np.empty_like(ells, dtype=float)\n",
      "for i,ell_max in enumerate(ells):\n",
      "    indices = np.array([[ell, mp, m] for ell in range(ell_max+1) for mp in range(-ell,ell+1) for m in range(-ell,ell+1)])\n",
      "    result = %timeit -o sp.wignerD(q.a, q.b, indices)\n",
      "    evals[i] = len(indices)\n",
      "    microseconds[i] = 1e6*result.best/len(indices)\n",
      "    print(\"With ell_max={0}, and {1} evaluations, each D component averages {2:.3} \u00b5s\".format(ell_max,evals[i],microseconds[i]))\n",
      "    sys.stdout.flush()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10000 loops, best of 3: 48.5 \u00b5s per loop\n",
        "With ell_max=1, and 10 evaluations, each D component averages 4.85 \u00b5s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10000 loops, best of 3: 151 \u00b5s per loop\n",
        "With ell_max=2, and 35 evaluations, each D component averages 4.32 \u00b5s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 351 \u00b5s per loop\n",
        "With ell_max=3, and 84 evaluations, each D component averages 4.18 \u00b5s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 675 \u00b5s per loop\n",
        "With ell_max=4, and 165 evaluations, each D component averages 4.09 \u00b5s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 1.18 ms per loop\n",
        "With ell_max=5, and 286 evaluations, each D component averages 4.13 \u00b5s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 1.85 ms per loop\n",
        "With ell_max=6, and 455 evaluations, each D component averages 4.07 \u00b5s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 2.6 ms per loop\n",
        "With ell_max=7, and 680 evaluations, each D component averages 3.82 \u00b5s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 3.73 ms per loop\n",
        "With ell_max=8, and 969 evaluations, each D component averages 3.85 \u00b5s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 5.12 ms per loop\n",
        "With ell_max=9, and 1330 evaluations, each D component averages 3.85 \u00b5s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 6.78 ms per loop\n",
        "With ell_max=10, and 1771 evaluations, each D component averages 3.83 \u00b5s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 8.8 ms per loop\n",
        "With ell_max=11, and 2300 evaluations, each D component averages 3.83 \u00b5s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 11.1 ms per loop\n",
        "With ell_max=12, and 2925 evaluations, each D component averages 3.79 \u00b5s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 14.4 ms per loop\n",
        "With ell_max=13, and 3654 evaluations, each D component averages 3.95 \u00b5s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 17.4 ms per loop\n",
        "With ell_max=14, and 4495 evaluations, each D component averages 3.86 \u00b5s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10 loops, best of 3: 20.3 ms per loop\n",
        "With ell_max=15, and 5456 evaluations, each D component averages 3.73 \u00b5s\n"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(evals, microseconds)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 22,
       "text": [
        "[<matplotlib.lines.Line2D at 0x102fb5690>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "svg": [
        "<?xml version=\"1.0\" encoding=\"utf-8\" standalone=\"no\"?>\n",
        "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
        "  \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
        "<!-- Created with matplotlib (http://matplotlib.org/) -->\n",
        "<svg height=\"641pt\" version=\"1.1\" viewBox=\"0 0 1065 641\" width=\"1065pt\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
        " <defs>\n",
        "  <style type=\"text/css\">\n",
        "*{stroke-linecap:butt;stroke-linejoin:round;}\n",
        "  </style>\n",
        " </defs>\n",
        " <g id=\"figure_1\">\n",
        "  <g id=\"patch_1\">\n",
        "   <path d=\"\n",
        "M0 641.477\n",
        "L1065.84 641.477\n",
        "L1065.84 0\n",
        "L0 0\n",
        "z\n",
        "\" style=\"fill:none;\"/>\n",
        "  </g>\n",
        "  <g id=\"axes_1\">\n",
        "   <g id=\"patch_2\">\n",
        "    <path d=\"\n",
        "M34.965 614.792\n",
        "L1039.37 614.792\n",
        "L1039.37 14.9425\n",
        "L34.965 14.9425\n",
        "z\n",
        "\" style=\"fill:#ffffff;\"/>\n",
        "   </g>\n",
        "   <g id=\"line2d_1\">\n",
        "    <path clip-path=\"url(#p15c525d802)\" d=\"\n",
        "M36.639 78.1022\n",
        "L40.824 307.721\n",
        "L49.0266 365.569\n",
        "L62.586 403.323\n",
        "L82.8414 387.741\n",
        "L111.132 414.652\n",
        "L148.797 519.113\n",
        "L197.176 506.779\n",
        "L257.607 508.652\n",
        "L331.43 517.005\n",
        "L419.985 518.04\n",
        "L524.61 531.527\n",
        "L646.645 466.78\n",
        "L787.428 501.811\n",
        "L948.299 559.788\" style=\"fill:none;stroke:#0000ff;stroke-linecap:square;\"/>\n",
        "   </g>\n",
        "   <g id=\"patch_3\">\n",
        "    <path d=\"\n",
        "M34.965 14.9425\n",
        "L1039.37 14.9425\" style=\"fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;\"/>\n",
        "   </g>\n",
        "   <g id=\"patch_4\">\n",
        "    <path d=\"\n",
        "M1039.37 614.792\n",
        "L1039.37 14.9425\" style=\"fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;\"/>\n",
        "   </g>\n",
        "   <g id=\"patch_5\">\n",
        "    <path d=\"\n",
        "M34.965 614.792\n",
        "L1039.37 614.792\" style=\"fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;\"/>\n",
        "   </g>\n",
        "   <g id=\"patch_6\">\n",
        "    <path d=\"\n",
        "M34.965 614.792\n",
        "L34.965 14.9425\" style=\"fill:none;stroke:#000000;stroke-linecap:square;stroke-linejoin:miter;\"/>\n",
        "   </g>\n",
        "   <g id=\"matplotlib.axis_1\">\n",
        "    <g id=\"xtick_1\">\n",
        "     <g id=\"line2d_2\">\n",
        "      <defs>\n",
        "       <path d=\"\n",
        "M0 0\n",
        "L0 -4\" id=\"m93b0483c22\" style=\"stroke:#000000;stroke-width:0.5;\"/>\n",
        "      </defs>\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"34.965\" xlink:href=\"#m93b0483c22\" y=\"614.7925\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"line2d_3\">\n",
        "      <defs>\n",
        "       <path d=\"\n",
        "M0 0\n",
        "L0 4\" id=\"m741efc42ff\" style=\"stroke:#000000;stroke-width:0.5;\"/>\n",
        "      </defs>\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"34.965\" xlink:href=\"#m741efc42ff\" y=\"14.9425\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"text_1\">\n",
        "      <!-- 0 -->\n",
        "      <defs>\n",
        "       <path d=\"\n",
        "M31.7812 66.4062\n",
        "Q24.1719 66.4062 20.3281 58.9062\n",
        "Q16.5 51.4219 16.5 36.375\n",
        "Q16.5 21.3906 20.3281 13.8906\n",
        "Q24.1719 6.39062 31.7812 6.39062\n",
        "Q39.4531 6.39062 43.2812 13.8906\n",
        "Q47.125 21.3906 47.125 36.375\n",
        "Q47.125 51.4219 43.2812 58.9062\n",
        "Q39.4531 66.4062 31.7812 66.4062\n",
        "M31.7812 74.2188\n",
        "Q44.0469 74.2188 50.5156 64.5156\n",
        "Q56.9844 54.8281 56.9844 36.375\n",
        "Q56.9844 17.9688 50.5156 8.26562\n",
        "Q44.0469 -1.42188 31.7812 -1.42188\n",
        "Q19.5312 -1.42188 13.0625 8.26562\n",
        "Q6.59375 17.9688 6.59375 36.375\n",
        "Q6.59375 54.8281 13.0625 64.5156\n",
        "Q19.5312 74.2188 31.7812 74.2188\" id=\"BitstreamVeraSans-Roman-30\"/>\n",
        "      </defs>\n",
        "      <g transform=\"translate(30.93375 630.95)scale(0.16 -0.16)\">\n",
        "       <use xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "    </g>\n",
        "    <g id=\"xtick_2\">\n",
        "     <g id=\"line2d_4\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"202.365\" xlink:href=\"#m93b0483c22\" y=\"614.7925\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"line2d_5\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"202.365\" xlink:href=\"#m741efc42ff\" y=\"14.9425\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"text_2\">\n",
        "      <!-- 1000 -->\n",
        "      <defs>\n",
        "       <path d=\"\n",
        "M12.4062 8.29688\n",
        "L28.5156 8.29688\n",
        "L28.5156 63.9219\n",
        "L10.9844 60.4062\n",
        "L10.9844 69.3906\n",
        "L28.4219 72.9062\n",
        "L38.2812 72.9062\n",
        "L38.2812 8.29688\n",
        "L54.3906 8.29688\n",
        "L54.3906 0\n",
        "L12.4062 0\n",
        "z\n",
        "\" id=\"BitstreamVeraSans-Roman-31\"/>\n",
        "      </defs>\n",
        "      <g transform=\"translate(183.415 630.95)scale(0.16 -0.16)\">\n",
        "       <use xlink:href=\"#BitstreamVeraSans-Roman-31\"/>\n",
        "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
        "       <use x=\"127.24609375\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
        "       <use x=\"190.869140625\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "    </g>\n",
        "    <g id=\"xtick_3\">\n",
        "     <g id=\"line2d_6\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"369.765\" xlink:href=\"#m93b0483c22\" y=\"614.7925\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"line2d_7\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"369.765\" xlink:href=\"#m741efc42ff\" y=\"14.9425\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"text_3\">\n",
        "      <!-- 2000 -->\n",
        "      <defs>\n",
        "       <path d=\"\n",
        "M19.1875 8.29688\n",
        "L53.6094 8.29688\n",
        "L53.6094 0\n",
        "L7.32812 0\n",
        "L7.32812 8.29688\n",
        "Q12.9375 14.1094 22.625 23.8906\n",
        "Q32.3281 33.6875 34.8125 36.5312\n",
        "Q39.5469 41.8438 41.4219 45.5312\n",
        "Q43.3125 49.2188 43.3125 52.7812\n",
        "Q43.3125 58.5938 39.2344 62.25\n",
        "Q35.1562 65.9219 28.6094 65.9219\n",
        "Q23.9688 65.9219 18.8125 64.3125\n",
        "Q13.6719 62.7031 7.8125 59.4219\n",
        "L7.8125 69.3906\n",
        "Q13.7656 71.7812 18.9375 73\n",
        "Q24.125 74.2188 28.4219 74.2188\n",
        "Q39.75 74.2188 46.4844 68.5469\n",
        "Q53.2188 62.8906 53.2188 53.4219\n",
        "Q53.2188 48.9219 51.5312 44.8906\n",
        "Q49.8594 40.875 45.4062 35.4062\n",
        "Q44.1875 33.9844 37.6406 27.2188\n",
        "Q31.1094 20.4531 19.1875 8.29688\" id=\"BitstreamVeraSans-Roman-32\"/>\n",
        "      </defs>\n",
        "      <g transform=\"translate(350.5225 630.95)scale(0.16 -0.16)\">\n",
        "       <use xlink:href=\"#BitstreamVeraSans-Roman-32\"/>\n",
        "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
        "       <use x=\"127.24609375\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
        "       <use x=\"190.869140625\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "    </g>\n",
        "    <g id=\"xtick_4\">\n",
        "     <g id=\"line2d_8\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"537.165\" xlink:href=\"#m93b0483c22\" y=\"614.7925\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"line2d_9\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"537.165\" xlink:href=\"#m741efc42ff\" y=\"14.9425\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"text_4\">\n",
        "      <!-- 3000 -->\n",
        "      <defs>\n",
        "       <path d=\"\n",
        "M40.5781 39.3125\n",
        "Q47.6562 37.7969 51.625 33\n",
        "Q55.6094 28.2188 55.6094 21.1875\n",
        "Q55.6094 10.4062 48.1875 4.48438\n",
        "Q40.7656 -1.42188 27.0938 -1.42188\n",
        "Q22.5156 -1.42188 17.6562 -0.515625\n",
        "Q12.7969 0.390625 7.625 2.20312\n",
        "L7.625 11.7188\n",
        "Q11.7188 9.32812 16.5938 8.10938\n",
        "Q21.4844 6.89062 26.8125 6.89062\n",
        "Q36.0781 6.89062 40.9375 10.5469\n",
        "Q45.7969 14.2031 45.7969 21.1875\n",
        "Q45.7969 27.6406 41.2812 31.2656\n",
        "Q36.7656 34.9062 28.7188 34.9062\n",
        "L20.2188 34.9062\n",
        "L20.2188 43.0156\n",
        "L29.1094 43.0156\n",
        "Q36.375 43.0156 40.2344 45.9219\n",
        "Q44.0938 48.8281 44.0938 54.2969\n",
        "Q44.0938 59.9062 40.1094 62.9062\n",
        "Q36.1406 65.9219 28.7188 65.9219\n",
        "Q24.6562 65.9219 20.0156 65.0312\n",
        "Q15.375 64.1562 9.8125 62.3125\n",
        "L9.8125 71.0938\n",
        "Q15.4375 72.6562 20.3438 73.4375\n",
        "Q25.25 74.2188 29.5938 74.2188\n",
        "Q40.8281 74.2188 47.3594 69.1094\n",
        "Q53.9062 64.0156 53.9062 55.3281\n",
        "Q53.9062 49.2656 50.4375 45.0938\n",
        "Q46.9688 40.9219 40.5781 39.3125\" id=\"BitstreamVeraSans-Roman-33\"/>\n",
        "      </defs>\n",
        "      <g transform=\"translate(517.94625 630.95)scale(0.16 -0.16)\">\n",
        "       <use xlink:href=\"#BitstreamVeraSans-Roman-33\"/>\n",
        "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
        "       <use x=\"127.24609375\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
        "       <use x=\"190.869140625\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "    </g>\n",
        "    <g id=\"xtick_5\">\n",
        "     <g id=\"line2d_10\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"704.565\" xlink:href=\"#m93b0483c22\" y=\"614.7925\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"line2d_11\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"704.565\" xlink:href=\"#m741efc42ff\" y=\"14.9425\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"text_5\">\n",
        "      <!-- 4000 -->\n",
        "      <defs>\n",
        "       <path d=\"\n",
        "M37.7969 64.3125\n",
        "L12.8906 25.3906\n",
        "L37.7969 25.3906\n",
        "z\n",
        "\n",
        "M35.2031 72.9062\n",
        "L47.6094 72.9062\n",
        "L47.6094 25.3906\n",
        "L58.0156 25.3906\n",
        "L58.0156 17.1875\n",
        "L47.6094 17.1875\n",
        "L47.6094 0\n",
        "L37.7969 0\n",
        "L37.7969 17.1875\n",
        "L4.89062 17.1875\n",
        "L4.89062 26.7031\n",
        "z\n",
        "\" id=\"BitstreamVeraSans-Roman-34\"/>\n",
        "      </defs>\n",
        "      <g transform=\"translate(685.1275 630.95)scale(0.16 -0.16)\">\n",
        "       <use xlink:href=\"#BitstreamVeraSans-Roman-34\"/>\n",
        "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
        "       <use x=\"127.24609375\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
        "       <use x=\"190.869140625\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "    </g>\n",
        "    <g id=\"xtick_6\">\n",
        "     <g id=\"line2d_12\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"871.965\" xlink:href=\"#m93b0483c22\" y=\"614.7925\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"line2d_13\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"871.965\" xlink:href=\"#m741efc42ff\" y=\"14.9425\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"text_6\">\n",
        "      <!-- 5000 -->\n",
        "      <defs>\n",
        "       <path d=\"\n",
        "M10.7969 72.9062\n",
        "L49.5156 72.9062\n",
        "L49.5156 64.5938\n",
        "L19.8281 64.5938\n",
        "L19.8281 46.7344\n",
        "Q21.9688 47.4688 24.1094 47.8281\n",
        "Q26.2656 48.1875 28.4219 48.1875\n",
        "Q40.625 48.1875 47.75 41.5\n",
        "Q54.8906 34.8125 54.8906 23.3906\n",
        "Q54.8906 11.625 47.5625 5.09375\n",
        "Q40.2344 -1.42188 26.9062 -1.42188\n",
        "Q22.3125 -1.42188 17.5469 -0.640625\n",
        "Q12.7969 0.140625 7.71875 1.70312\n",
        "L7.71875 11.625\n",
        "Q12.1094 9.23438 16.7969 8.0625\n",
        "Q21.4844 6.89062 26.7031 6.89062\n",
        "Q35.1562 6.89062 40.0781 11.3281\n",
        "Q45.0156 15.7656 45.0156 23.3906\n",
        "Q45.0156 31 40.0781 35.4375\n",
        "Q35.1562 39.8906 26.7031 39.8906\n",
        "Q22.75 39.8906 18.8125 39.0156\n",
        "Q14.8906 38.1406 10.7969 36.2812\n",
        "z\n",
        "\" id=\"BitstreamVeraSans-Roman-35\"/>\n",
        "      </defs>\n",
        "      <g transform=\"translate(852.75375 630.95)scale(0.16 -0.16)\">\n",
        "       <use xlink:href=\"#BitstreamVeraSans-Roman-35\"/>\n",
        "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
        "       <use x=\"127.24609375\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
        "       <use x=\"190.869140625\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "    </g>\n",
        "    <g id=\"xtick_7\">\n",
        "     <g id=\"line2d_14\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"1039.365\" xlink:href=\"#m93b0483c22\" y=\"614.7925\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"line2d_15\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"1039.365\" xlink:href=\"#m741efc42ff\" y=\"14.9425\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"text_7\">\n",
        "      <!-- 6000 -->\n",
        "      <defs>\n",
        "       <path d=\"\n",
        "M33.0156 40.375\n",
        "Q26.375 40.375 22.4844 35.8281\n",
        "Q18.6094 31.2969 18.6094 23.3906\n",
        "Q18.6094 15.5312 22.4844 10.9531\n",
        "Q26.375 6.39062 33.0156 6.39062\n",
        "Q39.6562 6.39062 43.5312 10.9531\n",
        "Q47.4062 15.5312 47.4062 23.3906\n",
        "Q47.4062 31.2969 43.5312 35.8281\n",
        "Q39.6562 40.375 33.0156 40.375\n",
        "M52.5938 71.2969\n",
        "L52.5938 62.3125\n",
        "Q48.875 64.0625 45.0938 64.9844\n",
        "Q41.3125 65.9219 37.5938 65.9219\n",
        "Q27.8281 65.9219 22.6719 59.3281\n",
        "Q17.5312 52.7344 16.7969 39.4062\n",
        "Q19.6719 43.6562 24.0156 45.9219\n",
        "Q28.375 48.1875 33.5938 48.1875\n",
        "Q44.5781 48.1875 50.9531 41.5156\n",
        "Q57.3281 34.8594 57.3281 23.3906\n",
        "Q57.3281 12.1562 50.6875 5.35938\n",
        "Q44.0469 -1.42188 33.0156 -1.42188\n",
        "Q20.3594 -1.42188 13.6719 8.26562\n",
        "Q6.98438 17.9688 6.98438 36.375\n",
        "Q6.98438 53.6562 15.1875 63.9375\n",
        "Q23.3906 74.2188 37.2031 74.2188\n",
        "Q40.9219 74.2188 44.7031 73.4844\n",
        "Q48.4844 72.75 52.5938 71.2969\" id=\"BitstreamVeraSans-Roman-36\"/>\n",
        "      </defs>\n",
        "      <g transform=\"translate(1020.095 630.95)scale(0.16 -0.16)\">\n",
        "       <use xlink:href=\"#BitstreamVeraSans-Roman-36\"/>\n",
        "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
        "       <use x=\"127.24609375\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
        "       <use x=\"190.869140625\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "    </g>\n",
        "   </g>\n",
        "   <g id=\"matplotlib.axis_2\">\n",
        "    <g id=\"ytick_1\">\n",
        "     <g id=\"line2d_16\">\n",
        "      <defs>\n",
        "       <path d=\"\n",
        "M0 0\n",
        "L4 0\" id=\"m728421d6d4\" style=\"stroke:#000000;stroke-width:0.5;\"/>\n",
        "      </defs>\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"34.965\" xlink:href=\"#m728421d6d4\" y=\"614.7925\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"line2d_17\">\n",
        "      <defs>\n",
        "       <path d=\"\n",
        "M0 0\n",
        "L-4 0\" id=\"mcb0005524f\" style=\"stroke:#000000;stroke-width:0.5;\"/>\n",
        "      </defs>\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"1039.365\" xlink:href=\"#mcb0005524f\" y=\"614.7925\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"text_8\">\n",
        "      <!-- 3.6 -->\n",
        "      <defs>\n",
        "       <path d=\"\n",
        "M10.6875 12.4062\n",
        "L21 12.4062\n",
        "L21 0\n",
        "L10.6875 0\n",
        "z\n",
        "\" id=\"BitstreamVeraSans-Roman-2e\"/>\n",
        "      </defs>\n",
        "      <g transform=\"translate(7.7475 619.2075)scale(0.16 -0.16)\">\n",
        "       <use xlink:href=\"#BitstreamVeraSans-Roman-33\"/>\n",
        "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-2e\"/>\n",
        "       <use x=\"95.41015625\" xlink:href=\"#BitstreamVeraSans-Roman-36\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "    </g>\n",
        "    <g id=\"ytick_2\">\n",
        "     <g id=\"line2d_18\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"34.965\" xlink:href=\"#m728421d6d4\" y=\"529.099642857\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"line2d_19\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"1039.365\" xlink:href=\"#mcb0005524f\" y=\"529.099642857\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"text_9\">\n",
        "      <!-- 3.8 -->\n",
        "      <defs>\n",
        "       <path d=\"\n",
        "M31.7812 34.625\n",
        "Q24.75 34.625 20.7188 30.8594\n",
        "Q16.7031 27.0938 16.7031 20.5156\n",
        "Q16.7031 13.9219 20.7188 10.1562\n",
        "Q24.75 6.39062 31.7812 6.39062\n",
        "Q38.8125 6.39062 42.8594 10.1719\n",
        "Q46.9219 13.9688 46.9219 20.5156\n",
        "Q46.9219 27.0938 42.8906 30.8594\n",
        "Q38.875 34.625 31.7812 34.625\n",
        "M21.9219 38.8125\n",
        "Q15.5781 40.375 12.0312 44.7188\n",
        "Q8.5 49.0781 8.5 55.3281\n",
        "Q8.5 64.0625 14.7188 69.1406\n",
        "Q20.9531 74.2188 31.7812 74.2188\n",
        "Q42.6719 74.2188 48.875 69.1406\n",
        "Q55.0781 64.0625 55.0781 55.3281\n",
        "Q55.0781 49.0781 51.5312 44.7188\n",
        "Q48 40.375 41.7031 38.8125\n",
        "Q48.8281 37.1562 52.7969 32.3125\n",
        "Q56.7812 27.4844 56.7812 20.5156\n",
        "Q56.7812 9.90625 50.3125 4.23438\n",
        "Q43.8438 -1.42188 31.7812 -1.42188\n",
        "Q19.7344 -1.42188 13.25 4.23438\n",
        "Q6.78125 9.90625 6.78125 20.5156\n",
        "Q6.78125 27.4844 10.7812 32.3125\n",
        "Q14.7969 37.1562 21.9219 38.8125\n",
        "M18.3125 54.3906\n",
        "Q18.3125 48.7344 21.8438 45.5625\n",
        "Q25.3906 42.3906 31.7812 42.3906\n",
        "Q38.1406 42.3906 41.7188 45.5625\n",
        "Q45.3125 48.7344 45.3125 54.3906\n",
        "Q45.3125 60.0625 41.7188 63.2344\n",
        "Q38.1406 66.4062 31.7812 66.4062\n",
        "Q25.3906 66.4062 21.8438 63.2344\n",
        "Q18.3125 60.0625 18.3125 54.3906\" id=\"BitstreamVeraSans-Roman-38\"/>\n",
        "      </defs>\n",
        "      <g transform=\"translate(7.835 533.514642857)scale(0.16 -0.16)\">\n",
        "       <use xlink:href=\"#BitstreamVeraSans-Roman-33\"/>\n",
        "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-2e\"/>\n",
        "       <use x=\"95.41015625\" xlink:href=\"#BitstreamVeraSans-Roman-38\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "    </g>\n",
        "    <g id=\"ytick_3\">\n",
        "     <g id=\"line2d_20\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"34.965\" xlink:href=\"#m728421d6d4\" y=\"443.406785714\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"line2d_21\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"1039.365\" xlink:href=\"#mcb0005524f\" y=\"443.406785714\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"text_10\">\n",
        "      <!-- 4.0 -->\n",
        "      <g transform=\"translate(7.365 447.821785714)scale(0.16 -0.16)\">\n",
        "       <use xlink:href=\"#BitstreamVeraSans-Roman-34\"/>\n",
        "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-2e\"/>\n",
        "       <use x=\"95.41015625\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "    </g>\n",
        "    <g id=\"ytick_4\">\n",
        "     <g id=\"line2d_22\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"34.965\" xlink:href=\"#m728421d6d4\" y=\"357.713928571\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"line2d_23\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"1039.365\" xlink:href=\"#mcb0005524f\" y=\"357.713928571\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"text_11\">\n",
        "      <!-- 4.2 -->\n",
        "      <g transform=\"translate(7.905 362.128928571)scale(0.16 -0.16)\">\n",
        "       <use xlink:href=\"#BitstreamVeraSans-Roman-34\"/>\n",
        "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-2e\"/>\n",
        "       <use x=\"95.41015625\" xlink:href=\"#BitstreamVeraSans-Roman-32\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "    </g>\n",
        "    <g id=\"ytick_5\">\n",
        "     <g id=\"line2d_24\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"34.965\" xlink:href=\"#m728421d6d4\" y=\"272.021071429\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"line2d_25\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"1039.365\" xlink:href=\"#mcb0005524f\" y=\"272.021071429\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"text_12\">\n",
        "      <!-- 4.4 -->\n",
        "      <g transform=\"translate(7.2 276.436071429)scale(0.16 -0.16)\">\n",
        "       <use xlink:href=\"#BitstreamVeraSans-Roman-34\"/>\n",
        "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-2e\"/>\n",
        "       <use x=\"95.41015625\" xlink:href=\"#BitstreamVeraSans-Roman-34\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "    </g>\n",
        "    <g id=\"ytick_6\">\n",
        "     <g id=\"line2d_26\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"34.965\" xlink:href=\"#m728421d6d4\" y=\"186.328214286\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"line2d_27\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"1039.365\" xlink:href=\"#mcb0005524f\" y=\"186.328214286\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"text_13\">\n",
        "      <!-- 4.6 -->\n",
        "      <g transform=\"translate(7.31 190.743214286)scale(0.16 -0.16)\">\n",
        "       <use xlink:href=\"#BitstreamVeraSans-Roman-34\"/>\n",
        "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-2e\"/>\n",
        "       <use x=\"95.41015625\" xlink:href=\"#BitstreamVeraSans-Roman-36\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "    </g>\n",
        "    <g id=\"ytick_7\">\n",
        "     <g id=\"line2d_28\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"34.965\" xlink:href=\"#m728421d6d4\" y=\"100.635357143\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"line2d_29\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"1039.365\" xlink:href=\"#mcb0005524f\" y=\"100.635357143\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"text_14\">\n",
        "      <!-- 4.8 -->\n",
        "      <g transform=\"translate(7.3975 105.050357143)scale(0.16 -0.16)\">\n",
        "       <use xlink:href=\"#BitstreamVeraSans-Roman-34\"/>\n",
        "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-2e\"/>\n",
        "       <use x=\"95.41015625\" xlink:href=\"#BitstreamVeraSans-Roman-38\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "    </g>\n",
        "    <g id=\"ytick_8\">\n",
        "     <g id=\"line2d_30\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"34.965\" xlink:href=\"#m728421d6d4\" y=\"14.9425\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"line2d_31\">\n",
        "      <g>\n",
        "       <use style=\"stroke:#000000;stroke-width:0.5;\" x=\"1039.365\" xlink:href=\"#mcb0005524f\" y=\"14.9425\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "     <g id=\"text_15\">\n",
        "      <!-- 5.0 -->\n",
        "      <g transform=\"translate(7.8175 19.3575)scale(0.16 -0.16)\">\n",
        "       <use xlink:href=\"#BitstreamVeraSans-Roman-35\"/>\n",
        "       <use x=\"63.623046875\" xlink:href=\"#BitstreamVeraSans-Roman-2e\"/>\n",
        "       <use x=\"95.41015625\" xlink:href=\"#BitstreamVeraSans-Roman-30\"/>\n",
        "      </g>\n",
        "     </g>\n",
        "    </g>\n",
        "   </g>\n",
        "  </g>\n",
        " </g>\n",
        " <defs>\n",
        "  <clipPath id=\"p15c525d802\">\n",
        "   <rect height=\"599.85\" width=\"1004.4\" x=\"34.965\" y=\"14.9425\"/>\n",
        "  </clipPath>\n",
        " </defs>\n",
        "</svg>\n"
       ],
       "text": [
        "<matplotlib.figure.Figure at 0x10ee0b150>"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "@njit #('void(complex128, complex128, int[:,:], complex128[:])')\n",
      "def _wignerD(Ra, Rb, indices, elements):\n",
      "    \"\"\"Main work function for computing Wigner D matrix elements\n",
      "\n",
      "    This is the core function that does all the work in the\n",
      "    computation, but it is strict about its input, and does not check\n",
      "    them for validity.\n",
      "\n",
      "    Input arguments\n",
      "    ===============\n",
      "    _wignerD(Ra, Rb, indices, elements)\n",
      "\n",
      "      * Ra, Rb are the complex components of the rotor\n",
      "      * indices is an array of integers [ell,mp,m]\n",
      "      * elements is an array of complex with length equal to the first\n",
      "        dimension of indices\n",
      "\n",
      "    The `elements` variable is needed because numba cannot create\n",
      "    arrays at the moment, but this is modified in place and returned.\n",
      "\n",
      "    \"\"\"\n",
      "    N = indices.shape[0]\n",
      "\n",
      "    # These constants are the recurring quantities in the computation\n",
      "    # of the matrix elements, so we calculate them here just once\n",
      "    absRa = abs(Ra)\n",
      "    absRb = abs(Rb)\n",
      "    absRRatioSquared = (absRb*absRb/(absRa*absRa) if absRa>epsilon else 0.0)\n",
      "    absRa_exp = (int(np.log(absRa)/_log2) if absRa>epsilon else min_exp)\n",
      "    absRb_exp = (int(np.log(absRb)/_log2) if absRb>epsilon else min_exp)\n",
      "\n",
      "    for i in range(N):\n",
      "        ell = indices[i,0]\n",
      "        mp = indices[i,1]\n",
      "        m = indices[i,2]\n",
      "\n",
      "        if(abs(mp)>ell or abs(m)>ell):\n",
      "            elements[i] = 0.0+0.0j\n",
      "\n",
      "        elif(absRa<=epsilon or 2*absRa_exp*(mp-m)<min_exp+mant_dig):\n",
      "            if mp!=m:\n",
      "                elements[i] = 0.0j\n",
      "            else:\n",
      "                if (ell+mp)%2==0:\n",
      "                    elements[i] = Rb**(2*m)\n",
      "                else:\n",
      "                    elements[i] = -Rb**(2*m)\n",
      "\n",
      "        elif(absRb<=epsilon or 2*absRb_exp*(mp-m)<min_exp+mant_dig):\n",
      "            if mp!=m:\n",
      "                elements[i] = 0.0j\n",
      "            else:\n",
      "                elements[i] = Ra**(2*m)\n",
      "\n",
      "        else:\n",
      "            rhoMin = max(0,mp-m)\n",
      "            rhoMax = min(ell+mp,ell-m)\n",
      "            if(absRa < 1.e-3):\n",
      "                # In this branch, we deal with NANs in certain cases when\n",
      "                # absRa is small by separating it out from the main sum,\n",
      "                Prefactor = coeff(ell, mp, m) * Ra**(m+mp) * Rb**(m-mp)\n",
      "                absRaSquared = absRa*absRa\n",
      "                absRbSquared = absRb*absRb\n",
      "                Sum = 0.0\n",
      "                for rho in range(rhoMax, rhoMin-1, -1):\n",
      "                    aTerm = absRaSquared**(ell-m-rho);\n",
      "                    if(aTerm != aTerm or aTerm<1.e-100):\n",
      "                        Sum *= absRbSquared\n",
      "                    else:\n",
      "                        if rho%2==0:\n",
      "                            Sum = ( binomial_coefficient(ell+mp,rho) * binomial_coefficient(ell-mp, ell-rho-m) * aTerm\n",
      "                                    + Sum * absRbSquared )\n",
      "                        else:\n",
      "                            Sum = ( -binomial_coefficient(ell+mp,rho) * binomial_coefficient(ell-mp, ell-rho-m) * aTerm\n",
      "                                    + Sum * absRbSquared )\n",
      "                elements[i] = Prefactor * Sum * absRbSquared**rhoMin\n",
      "            else:\n",
      "                Prefactor = coeff(ell, mp, m) * absRa**(2*ell-2*m) * Ra**(m+mp) * Rb**(m-mp)\n",
      "                Sum = 0.0\n",
      "                for rho in range(rhoMax, rhoMin-1, -1):\n",
      "                    if rho%2==0:\n",
      "                        Sum = (  binomial_coefficient(ell+mp,rho) * binomial_coefficient(ell-mp, ell-rho-m)\n",
      "                                 + Sum * absRRatioSquared )\n",
      "                    else:\n",
      "                        Sum = ( -binomial_coefficient(ell+mp,rho) * binomial_coefficient(ell-mp, ell-rho-m)\n",
      "                                 + Sum * absRRatioSquared )\n",
      "                elements[i] = Prefactor * Sum * absRRatioSquared**rhoMin\n",
      "\n",
      "#     return elements\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 77
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ell_max=8\n",
      "indices = np.array([[ell, mp, m] for ell in range(ell_max+1) for mp in range(-ell,ell+1) for m in range(-ell,ell+1)], dtype=int)\n",
      "elements = np.empty((indices.shape[0],), dtype=complex)\n",
      "elements_orig = np.copy(elements)\n",
      "_wignerD(q.a, q.b, indices, elements)\n",
      "#print(elements)\n",
      "print(np.allclose(elements, elements_orig))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "False\n"
       ]
      }
     ],
     "prompt_number": 78
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "indices"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 54,
       "text": [
        "array([[0, 0, 0]])"
       ]
      }
     ],
     "prompt_number": 54
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ell_max"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 55,
       "text": [
        "0"
       ]
      }
     ],
     "prompt_number": 55
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "elements, indices.shape[0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 43,
       "text": [
        "(array([  3.10503618e+231 -2.31584193e+77j]), 1)"
       ]
      }
     ],
     "prompt_number": 43
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "elements.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 46,
       "text": [
        "(1,)"
       ]
      }
     ],
     "prompt_number": 46
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "_wignerD.inspect_types()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "_wignerD (complex128, complex128, array(int32, 1d, A, nonconst), array(complex128, 1d, A, nonconst))\n",
        "--------------------------------------------------------------------------------\n",
        "# File: <ipython-input-11-9925f89de8f1>\n",
        "# --- LINE 5 --- \n",
        "\n",
        "@jit('void(c16, c16, i4[:], c16[:])')\n",
        "\n",
        "# --- LINE 6 --- \n",
        "\n",
        "def _wignerD(Ra, Rb, indices, elements):\n",
        "\n",
        "    # --- LINE 7 --- \n",
        "\n",
        "    \"\"\"Main work function for computing Wigner D matrix elements\n",
        "\n",
        "# --- LINE 8 --- \n",
        "\n",
        "\n",
        "\n",
        "    # --- LINE 9 --- \n",
        "\n",
        "    This is the core function that does all the work in the\n",
        "\n",
        "    # --- LINE 10 --- \n",
        "\n",
        "    computation, but it is strict about its input, and does not check\n",
        "\n",
        "    # --- LINE 11 --- \n",
        "\n",
        "    them for validity.\n",
        "\n",
        "# --- LINE 12 --- \n",
        "\n",
        "\n",
        "\n",
        "    # --- LINE 13 --- \n",
        "\n",
        "    Input arguments\n",
        "\n",
        "    # --- LINE 14 --- \n",
        "\n",
        "    ===============\n",
        "\n",
        "    # --- LINE 15 --- \n",
        "\n",
        "    _wignerD(Ra, Rb, indices, elements)\n",
        "\n",
        "# --- LINE 16 --- \n",
        "\n",
        "\n",
        "\n",
        "      # --- LINE 17 --- \n",
        "\n",
        "      * Ra, Rb are the complex components of the rotor\n",
        "\n",
        "      # --- LINE 18 --- \n",
        "\n",
        "      * indices is an array of integers [ell,mp,m]\n",
        "\n",
        "      # --- LINE 19 --- \n",
        "\n",
        "      * elements is an array of complex with length equal to the first\n",
        "\n",
        "        # --- LINE 20 --- \n",
        "\n",
        "        dimension of indices\n",
        "\n",
        "# --- LINE 21 --- \n",
        "\n",
        "\n",
        "\n",
        "    # --- LINE 22 --- \n",
        "\n",
        "    The `elements` variable is needed because numba cannot create\n",
        "\n",
        "    # --- LINE 23 --- \n",
        "\n",
        "    arrays at the moment, but this is modified in place and returned.\n",
        "\n",
        "# --- LINE 24 --- \n",
        "\n",
        "\n",
        "\n",
        "    # --- LINE 25 --- \n",
        "\n",
        "    \"\"\"\n",
        "\n",
        "# --- LINE 26 --- \n",
        "\n",
        "\n",
        "\n",
        "    # --- LINE 27 --- \n",
        "\n",
        "    # These constants are the recurring quantities in the computation\n",
        "\n",
        "    # --- LINE 28 --- \n",
        "\n",
        "    # of the matrix elements, so we calculate them here just once\n",
        "\n",
        "    # --- LINE 29 --- \n",
        "    # label 0\n",
        "    #   Ra.1 = Ra  :: pyobject\n",
        "    #   del Ra\n",
        "    #   Rb.1 = Rb  :: pyobject\n",
        "    #   del Rb\n",
        "    #   indices.1 = indices  :: pyobject\n",
        "    #   del indices\n",
        "    #   elements.1 = elements  :: pyobject\n",
        "    #   del elements\n",
        "    #   $0.1 = global(abs: <built-in function abs>)  :: pyobject\n",
        "    #   $0.3 = call $0.1(Ra.1, )  :: pyobject\n",
        "    #   del Ra.1\n",
        "    #   del $0.1\n",
        "    #   absRa = $0.3  :: pyobject\n",
        "    #   del $0.3\n",
        "\n",
        "    absRa = abs(Ra)\n",
        "\n",
        "    # --- LINE 30 --- \n",
        "    #   $0.4 = global(abs: <built-in function abs>)  :: pyobject\n",
        "    #   $0.6 = call $0.4(Rb.1, )  :: pyobject\n",
        "    #   del Rb.1\n",
        "    #   del $0.4\n",
        "    #   absRb = $0.6  :: pyobject\n",
        "    #   del $0.6\n",
        "\n",
        "    absRb = abs(Rb)\n",
        "\n",
        "    # --- LINE 31 --- \n",
        "    #   $0.8 = global(epsilon: 1e-14)  :: pyobject\n",
        "    #   $0.9 = absRa > $0.8  :: pyobject\n",
        "    #   del $0.8\n",
        "    #   branch $0.9, 36, 54\n",
        "    # label 36\n",
        "    #   del $0.9\n",
        "    #   $absRb36.1 = absRb  :: pyobject\n",
        "    #   $36.3 = $absRb36.1 * absRb  :: pyobject\n",
        "    #   del $absRb36.1\n",
        "    #   $absRa36.4 = absRa  :: pyobject\n",
        "    #   $36.6 = $absRa36.4 * absRa  :: pyobject\n",
        "    #   del $absRa36.4\n",
        "    #   $36.7 = $36.3 /? $36.6  :: pyobject\n",
        "    #   del $36.6\n",
        "    #   del $36.3\n",
        "    #   $phi57.1 = $36.7  :: pyobject\n",
        "    #   del $36.7\n",
        "    #   jump 57\n",
        "    # label 54\n",
        "    #   del $0.9\n",
        "    #   $const54.1 = const(float, 0.0)  :: pyobject\n",
        "    #   $phi57.1 = $const54.1  :: pyobject\n",
        "    #   del $const54.1\n",
        "    #   jump 57\n",
        "    # label 57\n",
        "    #   absRRatioSquared = $phi57.1  :: pyobject\n",
        "    #   del absRRatioSquared\n",
        "    #   del $phi57.1\n",
        "\n",
        "    absRRatioSquared = (absRb*absRb/(absRa*absRa) if absRa>epsilon else 0.0)\n",
        "\n",
        "    # --- LINE 32 --- \n",
        "    # label 97\n",
        "    #   del absRa\n",
        "    #   del $57.4\n",
        "    #   $97.1 = global(min_exp: -1021)  :: pyobject\n",
        "    #   $phi100.1 = $97.1  :: pyobject\n",
        "    #   del $97.1\n",
        "    #   jump 100\n",
        "    # label 100\n",
        "    #   absRa_exp = $phi100.1  :: pyobject\n",
        "    #   del absRa_exp\n",
        "    #   del $phi100.1\n",
        "    # label 72\n",
        "    #   del $57.4\n",
        "    #   $72.1 = global(int: <type 'int'>)  :: pyobject\n",
        "    #   $72.2 = global(np: <module 'numpy' from '/Users/boyle/.continuum/anaconda/lib/python2.7/site-packages/numpy/__init__.pyc'>)  :: pyobject\n",
        "    #   $72.3 = getattr(attr=log, value=$72.2)  :: pyobject\n",
        "    #   del $72.2\n",
        "    #   $72.5 = call $72.3(absRa, )  :: pyobject\n",
        "    #   del absRa\n",
        "    #   del $72.3\n",
        "    #   $72.6 = global(_log2: 0.69314718056)  :: pyobject\n",
        "    #   $72.7 = $72.5 /? $72.6  :: pyobject\n",
        "    #   del $72.6\n",
        "    #   del $72.5\n",
        "    #   $72.8 = call $72.1($72.7, )  :: pyobject\n",
        "    #   del $72.7\n",
        "    #   del $72.1\n",
        "    #   $phi100.1 = $72.8  :: pyobject\n",
        "    #   del $72.8\n",
        "    #   jump 100\n",
        "    #   $57.3 = global(epsilon: 1e-14)  :: pyobject\n",
        "    #   $57.4 = absRa > $57.3  :: pyobject\n",
        "    #   del $57.3\n",
        "    #   branch $57.4, 72, 97\n",
        "\n",
        "    absRa_exp = (int(np.log(absRa)/_log2) if absRa>epsilon else min_exp)\n",
        "\n",
        "    # --- LINE 33 --- \n",
        "    #   $100.3 = global(epsilon: 1e-14)  :: pyobject\n",
        "    #   $100.4 = absRb > $100.3  :: pyobject\n",
        "    #   del $100.3\n",
        "    #   branch $100.4, 115, 140\n",
        "    # label 140\n",
        "    #   del absRb\n",
        "    #   del $100.4\n",
        "    #   $140.1 = global(min_exp: -1021)  :: pyobject\n",
        "    #   $phi143.1 = $140.1  :: pyobject\n",
        "    #   del $140.1\n",
        "    #   jump 143\n",
        "    # label 143\n",
        "    #   absRb_exp = $phi143.1  :: pyobject\n",
        "    #   del absRb_exp\n",
        "    #   del $phi143.1\n",
        "    # label 115\n",
        "    #   del $100.4\n",
        "    #   $115.1 = global(int: <type 'int'>)  :: pyobject\n",
        "    #   $115.2 = global(np: <module 'numpy' from '/Users/boyle/.continuum/anaconda/lib/python2.7/site-packages/numpy/__init__.pyc'>)  :: pyobject\n",
        "    #   $115.3 = getattr(attr=log, value=$115.2)  :: pyobject\n",
        "    #   del $115.2\n",
        "    #   $115.5 = call $115.3(absRb, )  :: pyobject\n",
        "    #   del absRb\n",
        "    #   del $115.3\n",
        "    #   $115.6 = global(_log2: 0.69314718056)  :: pyobject\n",
        "    #   $115.7 = $115.5 /? $115.6  :: pyobject\n",
        "    #   del $115.6\n",
        "    #   del $115.5\n",
        "    #   $115.8 = call $115.1($115.7, )  :: pyobject\n",
        "    #   del $115.7\n",
        "    #   del $115.1\n",
        "    #   $phi143.1 = $115.8  :: pyobject\n",
        "    #   del $115.8\n",
        "    #   jump 143\n",
        "\n",
        "    absRb_exp = (int(np.log(absRb)/_log2) if absRb>epsilon else min_exp)\n",
        "\n",
        "# --- LINE 34 --- \n",
        "\n",
        "\n",
        "\n",
        "    # --- LINE 35 --- \n",
        "    # label 149.1\n",
        "    #   $const149.1.1 = const(LiftedLoop, LiftedLoop(<function _wignerD at 0x111f70aa0>))  :: pyobject\n",
        "    #   $149.1.4 = call $const149.1.1(indices.1, elements.1, )  :: pyobject\n",
        "    #   del indices.1\n",
        "    #   del elements.1\n",
        "    #   del $const149.1.1\n",
        "    #   del $149.1.4\n",
        "    #   jump 233\n",
        "    #   jump 149.1\n",
        "\n",
        "    for i,(ell,mp,m) in enumerate(indices):\n",
        "\n",
        "# --- LINE 36 --- \n",
        "\n",
        "\n",
        "\n",
        "        # --- LINE 37 --- \n",
        "\n",
        "        if(abs(mp)>ell or abs(m)>ell):\n",
        "\n",
        "            # --- LINE 38 --- \n",
        "            # label 233\n",
        "            #   $const233.1 = const(NoneType, None)  :: pyobject\n",
        "            #   $233.2 = cast(value=$const233.1)  :: pyobject\n",
        "            #   del $const233.1\n",
        "            #   return $233.2\n",
        "\n",
        "            elements[i] = 0.0+0.0j\n",
        "\n",
        "# The function contains lifted loops\n",
        "# Loop at line 35\n",
        "# Has 0 overloads\n",
        "\n",
        "================================================================================\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}